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Stability of the Singular Vortex and Associated Majorana Zero Modes in Trapped 
p-Wave Resonant Superfluids of Neutral Cold Atoms 
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The stability conditions for the singular vortex that accompanies Majorana zero modes at 
the core are investigated for p-wave resonant superfluids of atomic Fermi gases. Within the 
Ginzburg-Landau framework, we determine the stable conditions in the parameter space for 
the external rotation frequency and the harmonic trap frequency. There exists a narrow stable 
region in this parameter space for quasi-two-dimensional condensates. We also describe the 
detailed characterizations of the spatial structure of the order parameter in the chiral p-wave 
superfluids under rotation. 
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1. Introduction 

Chiral triplet superconductors and superfluids de- 
scribed symbolically by p + ip pairing, which breaks time 
reversal symmetry, are a typical topological order and 
possess many interesting and useful nontrivial physical 
properties if those really exist. As for superconductors, 
Sr2Ru04 is often referred to as a prime candidate chi- 
ral triplet superconductor. However, the precise pairing 
symmetry has not yet been identified and has been under 
strong discussion. A naive p x + ip y scenario is reexam- 
ined from various aspects. 1 - 1 Therefore, at this moment, 
there is no concrete material that firmly exhibits chiral 
triplet pairing symmetry in a superconductor. 2 ' 

As for superfluids, the 3 He-A phase is definitely in the 
chiral pairing state, which has been a testing ground for 
checking ideas associated with the topological order. 3 ) 
However, because of its inherent strong interacting na- 
ture from the theoretical side and the limited accessibility 
from the experimental side, it is a difficult and challeng- 
ing task to explore supernuid 3 He. Needless to say, we 
wish to find more clean systems of chiral superfluids that 
are easily accessible theoretically and experimentally. 

In this sense, the use of neutral atomic Fermi gases is 
another option to pursue this possibility. By using mag- 
netic Feshbach resonance, a p-wave interaction channel 
can be enhanced, leading to p-wave paired superfluids. In 
fact, there have been vigorous activities toward realizing 
p-wave resonant superfluids in 6 Li 4 ) and 40 K 5 ) recently. 
The advantages of p-wave resonant superfluids obtained 
using atomic Fermi gases are (1) the tunability of the 
attractive interaction strength in a p-wave channel from 
strong to weak cases, (2) the controllability of the sys- 
tem dimension by changing a trap potential from three 
dimensions (3D) to two dimensions (2D), even to one 
dimension (ID), and (3) ease of manipulations, such as 
rotation, and excitations of collective modes by external 
perturbations. 

Chiral p-wave superfluids have nontrivial Caroli-de 
Gennes-Matricon (CdGM) states 6 ) at a vortex core. It 
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is possible for the CdGM states to have exact zero en- 
ergy states in the BCS regime. 7 ~ 9 > Since the zero energy 
states can be described by the ID Majorana equation, 10 ) 
the zero energy states are Majorana zero modes. A re- 
markable feature of Majorana zero modes is that their 
creation operator is expressed by the self-Hermitian oper- 
ator 7 = 7^ and their host vortices obey the non-Abelian 
statistics. 11 ' 12 ) These nontrivial properties can be uti- 
lized for a topological quantum computer. 13 ) 

A pair of conventional fermionic creation and an- 
nihilation operators span a 2D Hilbert space because 
their square vanishes. This is not true for the Majo- 
rana operators. Thus, to avoid the problem when 2Nq 
vortices are present, we can construct "conventional" 
complex fermionic creation and annihilation operators, 
ip} = (li - ijN„+i)/2 and tpi = (7* + ij No+i )/2, respec- 
tively, from a pair of vortices, where the normalization 
is chosen as jf = 1 and subscript i denotes the index 
of vortices. 14 ^ These operators satisfy ipf — = and 
thus span a 2D subspace of degenerate ground states as- 
sociated with these operators. Strictly speaking, a topo- 
logical vortex qubit is defined through two pairs of vor- 
tices because the fermionic number is fixed to be even or 
odd. 15 ) 

A necessary condition for Majorana zero modes to ex- 
ist at a vortex core is as follows: (A) the vortex possesses 
an odd winding number and phase singularity in spinless 
chiral p-wave superfluids 9 ' 16 ) or (B) the vortex possesses 
one-half winding number, that is, the vortex is a half- 
quantum vortex in spinful chiral p-wave superfluids. 11 ) 

Here, under those circumstances, we investigate exper- 
imental requirements to satisfy the condition (A) in neu- 
tral atomic Fermi gases. Since the population of atoms 
in each hyperfine spin state is controllable by an RF 
field, the spinless p-wave superfluids could be achieved. In 
quasi-2D confinement that prohibits the Cooper pairs of 
the z-direction, p x ± ip y paired states are realized unless 
the magnetic field for the Feshbach resonance is parallel 
to the quasi-2D system. 17, 1S ) A requirement of the quasi- 
2D system is R z < £, namely, hio z /(k B T c ) > T C /T F , 
where R z is the radius of atomic gases toward the z- 
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direction, £ is the coherence length, co z is the trap fre- 
quency of the z-direction, T c is the transition tempera- 
ture, and Tf is the Fermi temperature. This requirement 
may be satisfied by the optical lattice potential of a pe- 
riod ~ 3 fim. 19 ' The aim of this study is to determine 
the stability condition for the singular vortex, which pos- 
sesses Majorana zero energy modes at the vortex core, 
under the quasi-2D confinement that allows the chiral 
p x ± ip y paired state produced by the magnetic Fesh- 
bach resonance. It has been proposed that the vortices 
with Majorana zero modes in the p-wave resonant atomic 
Fermi gases are utilized for a topological quantum com- 
puter already including the optical methods of operat- 
ing and reading out states of qubits. 8,15 ^ However, the 
strict stable conditions of the vortices with Majorana 
zero modes have not yet been discussed. 

In §2, we derive the Ginzburg-Landau (GL) functional 
form within the weak coupling approximation in order 
to describe the p x ± ip y paired state under rotation, and 
classify the possible vortex types in a 2D trapped system. 
We construct the phase diagram for the stable singular 
vortex, which allows Majorana zero modes in the core, 
in the parameter space, that is, rotation frequency vs 
trap frequency lj in §3. The final section is devoted to 
conclusions. The closely related vortex stability calcula- 
tions for 3D systems are found in refs. 20 and 21. 

2. Formulation 

2. 1 Ginzburg-Landau functional 

The order parameter (OP) of the spinless p-wave su- 
perfluidity in the quasi-2D trap potential can be ex- 
panded on the basis of rectangular coordinates with ex- 
pansion coefficients Ai (i = x,y): 



/ha 



-tu 2 r 2 A*A t 
5 



(6) 



A(P) = A X P X + Aypy, 



(1) 



where p is the unit vector in momentum space. We as- 
sume that the momentum toward the z-direction is pro- 
hibited by the quasi-2D trap potential. 

Here, we employ the GL framework for the p-wave su- 
perfluidity. 3 ) This framework is valid in the vicinity of 
the transition temperature, T c — T <C T c , and applicable 
to cold atomic Fermi gases in a harmonic trap potential 
with hw -C ksT c , 22 ^ where T c and to are the transition 
temperature and the trap frequency, respectively. The 
GL free energy is obtained by integrating the GL free 
energy density in a 2D system: 



F = J d 2 r(j 'bulk + ,/grad + /cent + /harm 



(2) 



In the weak coupling limit, the condensation energy, gra- 
dient energy, centrifugal energy, and trap potential en- 
ergy are respectively given by 

1 * 1 

/bulk = — (1 — t) A i Ai + —A i AiAjAj + —A i A i AjAj, 

(3) 

/grad = f[(3i^) Wi) + (d t A 3 y(d 3 AA 

+ (d i A i y(d j A j )}, (4) 
f cent = -ln 2 r 2 (A*A i + 2\A e \ 2 ), (5) 



where t = T/T c , d = V — ifl x r with the rota- 
tion frequency f2 around the z-axis, H | z, and Ag = 
— A x s'm9 + A y cos 9 in cylindrical coordinates. The re- 
peated indices i and j imply summations over x and y. 

The GL free energy density functional, eqs. (3) - 
(6), are expressed with dimensionlcss units. The units 
of OP and length are the OP amplitude of the chiral 
state in a bulk, A = v /107r 2 /7C(3)fc s T c , and the GL 
coherence length £ = / A%n 2 (hv F /k B T c ) with 

zero temperature, respectively, where vf is the Fermi 
velocity. The effect of harmonic trap potential on con- 
densates depends on the dimensionless pairing interac- 
tion r that determines the transition temperature T c = 
CEp exp(— 1/r), where Ep is the Fermi energy and C 
is a numerical coefficient of order unity. 22 ) We use the 
unit hcu = v /1447r 2 /35C(3) v /2r/(l + 2Y)k B T c for the 
trap frequency ftw. The p-wave triplet paring interaction 
is r w A 2 /13 with a gaseous parameter A = 2\a\pp /irh, 
where a is the scattering length and pp is the Fermi mo- 
mentum. 23 ) Rotation frequency is normalized by the crit- 
ical rotation frequency f2 c , where condensates fly away 
from the trap potential by centrifugal force. We sum up 
the centrifugal energy and the trap potential energy as 



/cent ~t~ /ha 



3, 2 



Q 2 )r 2 \A r \ 2 + Uco 2 -3n 2 )r 2 \A s \ 2 
5 



with cylindrical coordinates. Therefore, the critical rota- 
tion frequency is given by Q c — uj/^/3. 20 ' 21 ) 

We estimate the approximate profile of the conden- 
sates at rest within the Thomas-Fermi (TF) approxima- 
tion. In addition, we assume that the condensate is in 
the chiral state, which is stable in a bulk. The profile of 
OP amplitude is 



|A c (r)| = ^(l-i)- 
and the radius of condensates is 



(7) 



(8) 



Here, we introduce the Z-vector pointing in the direc- 
tion of the orbital angular momentum of the Cooper 
pair. The Z-vector has only the z-componcnt when the 
momentum toward the z-direction is suppressed. The 
z-component of the /-vector is expressed with expan- 
sion coefficients on the basis of the component of the 
orbital angular momentum along the z-axis: A(p) = 
A + p + + A_p_, where the bases and coefficients of the 
plus- minus chiral states are p± = T (p x ± ip y ) /\/2 and 
A± = =F (A x =F iA y ) /\/2, respectively. By using the co- 
efficients, 



\A+\ 2 - \A_ 



(9) 



where |A| 2 = \A + \ 2 + \A-\ 2 is the squared amplitude of 
the OP. If the condensates are purely in the chiral state, 
l z has the unit value. l z decreases when the polar state is 
mixed through gradient and centrifugal energies. Finally, 
l z becomes zero, provided that the condensates are in the 
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pure polar state. 

We have identified stationary solutions by numerically 
solving the variational equation Sf(r)/8Ai(r) = 0, where 
/(r) is the GL free energy density functional, namely, the 
integrand of eq. (2) . In this study, we have considered the 
temperature t = 0.9. 

2.2 Possible types of axisymmetric structures and vor- 
tices at low rotation 

Since we assume the axisymmetric trap potential, ax- 
isymmetric structures of the condensates are stable when 
there is no vortex and only a singular vortex exists in 
the system. Here, we define the winding numbers of 
the plus and minus components as w+ and respec- 
tively, where A± = \A±\e l( > w±e+a± ' > with a phase con- 
stant a± . Axisymmetric structures of the OP consisting 
of the plus-minus components must have the combina- 
tions of the winding number (ui+, ui_) — (n — 1, n + 1), 
where n is an integer. 24 ' 25 ) The axisymmetric condi- 
tion of what is stricter than that of multicomponent 
spinor Bose-Einstein condensates (BECs). 26 ) Generally, 
the small winding number state is energetically favorable. 

For the nonvortex structures, the possible combina- 
tions of the winding number are (u>+, u>_) = (—2,0) and 
(0,2). These structures are degenerate at rest. However, 
under positive rotation or counterclockwise rotation, the 
structure with (—2,0) is stable because the minus chi- 
ral component without phase winding is advantageous 
under positive rotation. 27 ' For the structures with a sin- 
gular vortex, the possible combinations of the winding 
number are (w + ,w-) = (—3,-1), (—1,1), and (1,3). 
Among these combinations, the structure with (—1, 1) is 
the most stable under low positive rotation, because the 
main component, which is the minus chiral component, 
has a winding number of 1. 

Therefore, it is sufficient to consider the combination 
of winding number (w + , w_) — (—2, 0) for the nonvortex 
structure, and (w + , w_) = (—1, 1) for the structure with 
a singular vortex. 

3. Stability Region of the Singular Vortex 

We show the phase diagram of the OP structures for 
trap frequency lo vs rotation frequency Q in Fig. 1. We 
consider the three cases for the trap frequencies hw = 2.5, 
4.0, and 8.0 x 10" V 1447r 7 35 C(3)\/2r/(l + 2T)k B T c . 
The condition of applying the GL framework huj <C ksT c 
is satisfied because of 2r/(l + 2r) < 1. The calculated 
points are denoted by solid circles and crosses in Fig. 1. 
The region N is the nonvortex structure, S is the singular 
vortex structure where Majorana zero modes exist, and 
M is the multiple vortex region. Details of each structure 
are explained in the following subsections. 

3.1 Rotation dependence 

In this subsection, we show the representative struc- 
tures with the trap frequency 2.5 x 10~ 3 woi where the 
TF radius of the condensate Ra is about 160£ at rest. 
In this case, various vortex structures are stabilized by 
rotation. 

We show the OP amplitude at rest in Figs. 2(a) and 
2(b). The OP amplitude at the center is slightly sup- 
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Fig. 1. (Color online) Phase diagram of the OP structures for 
trap frequency (uj) vs rotation frequency (f2). The TF radii of 
the condensate (Ra) at rest are also noted. The region N is the 
nonvortex structure, S is the singular vortex structure where Ma- 
jorana zero modes exist, and M is the multiple vortex structure. 
Solid circles and crosses show the structures with and without 
Majorana zero modes, respectively. 

pressed from that of the chiral state in a bulk A(i) = 
AoVl — t. The OP amplitude is gradually decreased to- 
ward the outside by the trap potential, according to eq. 
(7) except near the edge. Since the condensation energy 
is small at the edge, the gradient energy makes a large 
contribution to the structures. The plus chiral compo- 
nent with the winding number w + = — 2 is induced 
at the edge because the gradient energy terms in eq. (4), 
such as A+(d 2 A_/dr 2 ), inevitably induce another com- 
ponent at the place where the main component A- with 
the winding number w_ = is spatially varied. The fea- 
ture is absent in spinor BECs because the spatial differ- 
ential operators act on real space, while the OP exists in 
spin space. 28,29 ) 

This is also clear from the l z values shown in Figs. 
2(c) and 2(d). In almost all regions, l z = —1, namely, 
the condensate is in the minus chiral state. However, the 
/-vectors shorten toward the outside of the condensate at 
the edge. The condensate is finally in the polar state. 

A singular vortex enters the center of the condensate 
under rotation. We show the OP amplitude and the z- 
component of the /-vector in Fig. 3 at f2 = 0.35fi c . Since 
the vortex enters the condensate of the spinless chiral 
state and has a combination of odd winding number 
(u> + ,u;_) = (—1,1), the vortex is accompanied by Ma- 
jorana zero modes. The condensate is in the minus chi- 
ral state in most of the regions; however, the plus chiral 
component is induced around the vortex core and at the 
edge. Note that the condensate at the edge is in the plus 
chiral state because of the effect of rotation through the 
gradient energy terms in eq. (4), such as VlrA* + (dA_/dr). 

On increasing the rotation frequency fi, vortices enter 
the condensate one by one from outside. We show the 
OP amplitude in Fig. 4(a) and l z in Fig. 4(b) where the 
condensate contains three vortices at SI = 0.60S1 C . All the 
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Fig. 2. (Color online) In the region N of Fig. 1 at rest: (a) Contour 
map of the OP amplitude, (b) cross sections of the OP amplitude, 
the plus chiral component, and TF estimate from eq. (7), (c) 
contour map of l z where the region A/Ao > 0.01 is shown, and 
(d) cross section of l z . 
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Fig. 3. (Color online) In the region S of Fig. 1 where a singular 
vortex is accompanied by Majorana zero modes at O = 0.35f2 c : 
(a) Contour map of the OP amplitude, (b) cross sections of the 
OP amplitude (solid line) and the plus chiral component (dot- 
dash line), (c) contour map of l z where the region A/Ao > 0.01 
is shown, and (d) cross section of l z . 



vortices are not singular vortices in the system. The cores 
of the plus and minus chiral components are split and 
covered by the other component. The centrifugal energy, 
eq. (5), increases the condensate extension and brings the 
condensate at the edge into the polar state with the mo- 
mentum toward the ^-direction. Half-quantum vortices 
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Fig. 4. (Color online) In the region M of Fig. 1 at Q = 0.60f2 c : 
(a) Contour map of the OP amplitude and (b) contour map of 
l z where the region A/Aq > 0.01 is shown. 



with l z = 1 in the polar state effectively absorb the an- 
gular momentum by rotation. 

Here, we introduce a pseudo-spin representation so 
that differences between the singular vortex and the split- 
ting vortex can easily be seen. Since the OP has only two 
components, we can express the OP on the basis of the 
vector 

! Aa 



1 



A. 



(10) 



The pseudo-spin is defined using the vectors and Pauli 
matrices with the total spin 1/2. 



S x 



S a 



l 



|A| 2 
1 



Re(A* + A^ 



lm(A* + A-) 



(11) 



S -Vcrn 1 |A+ ' 2 



A_ 



The z-component of the pseudo-spin S z is one- half of l z . 

We define the topological charge N = u>_ — u>+. Ax- 
isymmetric vortices have the topological charge N = 2 
and the half-quantum vortices with the combinations of 
the winding number (w+, iu_) = (—1,0) and (1,0) have 
N = 1. The pseudo-spin can be rewritten with the topo- 
logical charge: 



\A+A- 
|A|2 



• cos(N9 - 



_1 L4 + | 2 -|, 
z ~2 IAP 



■a'), 
a'), 



(12) 



where the relative phase a' = a- — a+. The pseudo-spin 
vector (S x , S v ) rotates N times around a vortex with the 
topological charge ./V. 30 ) 

The axisymmetric singular vortex at the center of the 
condensate at = 0.35£l c has the topological charge 
N = 2 (Figs. 5(a) and 5(b)). Similarly, the vortices at 
£1 = 0.60£l c have the topological charge TV — 2 around a 
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Fig. 5. (Color online) Pseudo-spin vectors (S x , S y ) with the con- 
tour of the Sz component (a) around a singular vortex at CI = 
0.35Q C (Fig. 3) and (c) around a splitting vortex at f! = 0S0Q c 
(Fig. 4). Cross sections of the OP amplitude (solid line), the 
plus chiral component (dot-dash line), and minus chiral compo- 
nent (dashed line) (b) across a singular vortex and (d) across a 
splitting vortex. 

path far away from them (Figs. 5(c) and 5(d)). However, 
they are split into two half-quantum vortices with N = 1 
because of the OP amplitude around the off-centered vor- 
tices without axisymmetry by a harmonic trap potential. 
The half-quantum vortices on the y-axis with the com- 
binations of the winding number (—1,0) and (0,1) are 
located on (0,42£o) and (0, 50£o), respectively, in Fig. 
5(d). The distance between two half-quantum vortices is 
approximately 8£o ~ 3£(£), where £(f) = £o(l — t)~ x / 2 is 
the temperature-dependent coherence length. 

Since the off-centered vortex has a finite OP every- 
where, Majorana zero modes, which are zero energy exci- 
tations, may not exist. Note, however, that the vortex sit- 
uated at a few times of coherence length is similar to the 
singular vortex with Majorana zero modes. The quasi- 
particle states of the off-centered vortex are unknown at 
present. We must solve the Bogoliubov-de Gennes equa- 
tions in order to understand the quasi-particle states. 
This is our future problem. 

3.2 Trapping dependence 

The number of vortices is proportional to the rota- 
tion frequency and the system size, or the cross section 
that is perpendicular to the rotational axis. The radius of 
the condensate is inversely proportional to the trap fre- 
quency and is increased by the centrifugal energy under 
rotation. Since the effect of the centrifugal energy dif- 
fers among orbital states, the region of the polar state 
appears predominantly at the edge of the condensate 
under high rotation. Therefore, the number of vortices 
in the chiral state situated in the central region does 
not increase with increasing rotation because the angu- 
lar momentum is absorbed by half-quantum vortices in 
the polar state at the surrounding area, which prevent 
the vortex from entering the central region. 



The vortex does not enter the chiral state up to high 
rotation in the case of co = 8.0 x 10 _3 wo and only the 
singular vortex enters under high rotation in the case of 
u) = 4.0 x 10~ 3 wo- Namely, we will have to employ a low 
trap frequency in an experiment to observe Majorana 
zero modes. In the previous studies under a 3D pancake- 
shape trap potential, 20,21 ) the major radius of the con- 
densate was approximately 30£o, which corresponds to 
high-trap- frequency cases in the present study; therefore, 
the polar core vortex with a combination of winding num- 
ber (w+, wo, W-) = (—1,0,1), which corresponds to the 
singular vortex (w + , w—) = (—1, 1) in this paper, did not 
enter the condensate, where Wo is the winding number 
of the OP z-component, which is absent in the present 
study. 

4. Conclusions 

We have investigated the stability of the singular vor- 
tex that possesses Majorana zero modes at the core in 
the p-wave superfluidity of the atomic Fermi gases by 
the GL framework. As a result, in the systems confined 
in the quasi-2D harmonic trap potential with a low fre- 
quency, the singular vortex accompanied by Majorana 
zero modes enters the center of the condensates under 
rotation. Under high rotation, the multiple vortex state 
is stabilized. The off-centered vortices in the multiple vor- 
tex state are split into two half-quantum vortices and are 
not accompanied by Majorana zero modes; the quasi- 
particle states of the off-centered vortex are unknown at 
present. Note that a centered singular vortex is accompa- 
nied by Majorana zero modes even though some vortices 
exist in the chiral state. 

Within the present quasi-2D harmonic trap potential, 
there is no stable parameter region in the (w, fi) plane for 
plural singular vortices, which is necessary for a topologi- 
cal quantum computer. We suggest the use of a quasi-2D 
square well potential for confinement, which is also used 
in the superfluid 3 He confined in parallel plates. 31 ' 32 ) 
This could accommodate a pair of singular vortices or 
more. 
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